set.seed(0401)
options (warn = -1)
# packbirth loading
library(pacman)
p_load("here",
       "haven",
       "car",
       "tidyverse",
       "ggplot2",
       "bruceR",
       "tinytable",
       "modelsummary",
       "gtsummary",
       "scales",
       "dotwhisker",
       "broom",
       "brms",
       "rstanarm",
       "flextable",
       "officer",
       "broom.mixed")

load("fake.Rdata")

data_raw$age <- 2024 - data_raw$birthyear
data_raw$age_cate <- Recode(data_raw$age,"20:29 = 0;30:39 = 1;40:65 = 2")
data_raw$age_cate <- factor(data_raw$age_cate,
                            levels = c(0,1,2))
data_raw$socialLevel_cate <- Recode(data_raw$socialLevel,"1:5 = 1;6:10 = 0")
data_raw$partymember <- Recode(data_raw$partymember,"2 = 0")
data_raw$race <- Recode(data_raw$race,"2 = 0")
data_raw$nationlism <- data_raw %>% 
  dplyr::select(B1__1:B1__3) %>% 
  rowMeans(na.rm = TRUE)
data_raw$edu_cate <- Recode(data_raw$edu,"4:5 = 4")
data_raw$edu_cate <- factor(data_raw$edu_cate,
                            levels = c(1,2,3,4),
                            labels = c("Primary","Junior","High School","College"))
data_raw$edu_order <- factor(data_raw$edu,
                            levels = c(1,2,3,4,5),
                            labels = c("Primary","Junior","High School","College","Postgraduate"))
data_raw$edu_dummy <- Recode(data_raw$edu,"1:3 = 0;4:5 = 1")
data_raw$broad_dummy <- Recode(data_raw$broad,"1:3 = 1")
data_raw$exposure_mean_dummy <- if_else(data_raw$exposure >= mean(data_raw$exposure),1,0)
data_raw$exposure_median_dummy <- if_else(data_raw$exposure >= median(data_raw$exposure),1,0)
data_raw$exposure_scale <- scales::rescale(data_raw$exposure,to = c(0,1))



data_raw$E1ar <- if_else(data_raw$E1a == 3,1,0)
data_raw$E1br <- if_else(data_raw$E1b == 1,1,0)
data_raw$E2ar <- if_else(data_raw$E2a == 3,1,0)
data_raw$E2cr <- if_else(data_raw$E2c == 1,1,0)
data_raw$F1ar <- if_else(data_raw$F1a == 1,1,0)
data_raw$F1br <- if_else(data_raw$F1b == 2,1,0)
data_raw$F2ar <- if_else(data_raw$F2a == 1,1,0)
data_raw$F2br <- if_else(data_raw$F2b == 2,1,0)
data_raw <- data_raw %>%
  mutate(know_pol = rowSums(select(., E1ar:F2br), na.rm = TRUE))

nationalist_news <- c("G1", "G2", "G6", "G7", "G11", "G12", "G14", "G16", "G17", "G18", "G19")
internationalist_news <- c("G3","G4","G5","G8","G9","G10","G13","G15","G20")
true_news <- c("G1","G4","G5","G7","G10","G11","G14","G15","G19","G20")
fake_news <- c("G2","G3","G6","G8","G9","G12","G13","G16","G17","G18")
nationalist_true_news <- intersect(nationalist_news, true_news)
nationalist_fake_news <- setdiff(nationalist_news, nationalist_true_news)

internationalist_true_news <- intersect(internationalist_news, true_news)
internationalist_fake_news <- setdiff(internationalist_news, internationalist_true_news)

data_raw <- data_raw %>% 
  mutate(across(G1:G20, ~ ifelse(. == 2, 0, .)),
         accuracy_score = (
           (G1 == 1) + (G4 == 1) + (G5 == 1) + (G7 == 1) + (G11 == 1) +
             (G14 == 1) + (G15 == 1) + (G19 == 1) + (G20 == 1) +(G10 == 1) +
             (G2 == 0) + (G3 == 0) + (G6 == 0) + (G8 == 0) + (G9 == 0) + 
             (G12 == 0) + (G13 == 0) + (G16 == 0) + 
             (G17 == 0) + (G18 == 0)
    ),
    accuracy_rate = accuracy_score/20
  ) %>% 
  mutate(
    nationalist_fake_accuracy = rowSums(across(all_of(nationalist_fake_news)) == 0) / length(nationalist_fake_news),
    internationalist_fake_accuracy = rowSums(across(all_of(internationalist_fake_news)) == 0) / length(internationalist_fake_news))

data_raw_long <- data_raw %>% 
  pivot_longer(
    cols = G1:G20,           
    names_to = "news_id",   
    values_to = "news_judgement" 
  ) %>% 
  mutate(
    news_type = case_when(
    news_id %in% c("G1","G4","G5","G7","G10","G11","G14","G15","G19","G20") ~ "1",
    TRUE ~ "0"),
    news_correct = if_else(news_judgement == news_type,1,0),
    news_category = case_when(
      news_id %in% nationalist_true_news ~ "Nationalist & True",
      news_id %in% nationalist_fake_news ~ "Nationalist & Fake",
      news_id %in% internationalist_true_news ~ "Internationalist & True",
      news_id %in% internationalist_fake_news ~ "Internationalist & Fake",
      TRUE ~ "Uncategorized"  # 以防有未分类的题目
    )
  )
data_raw_long$news_category <- factor(data_raw_long$news_category,levels = c("Nationalist & True","Nationalist & Fake","Internationalist & True","Internationalist & Fake"))



df_correct_rates <- data_raw_long %>%
  group_by(news_category) %>%
  summarise(
    correct_rate = mean(news_correct, na.rm = TRUE),
    .groups = "drop"
  ) %>% 
  mutate(overall_accuracy = mean(data_raw_long$news_correct, na.rm = TRUE))

p1 <- ggplot(df_correct_rates, aes(x = news_category, y = correct_rate)) +
  geom_bar(position = position_dodge(),stat = "identity", width = 0.5, fill = "darkgray") +
  scale_y_continuous("", limits = c(0, 1), label = percent) + 
  geom_hline(yintercept = df_correct_rates$overall_accuracy[1], 
             linetype = "dashed", color = "black", linewidth = .5) +
  labs(title = "Accuracy in News Judgement",
       x = "News Category", 
       y = "Correct Response Rate") +
  scale_x_discrete(labels = c(
    "Nationalist & True" = "China-praising(True)",
    "Nationalist & Fake" = "China-praising(Fake)",
    "Internationalist & True" = "China-disparaging(True)",
    "Internationalist & Fake" = "China-disparaging(Fake)"
  )) +
  theme_bw(base_size = 5) + 
  theme(
    panel.background = element_rect(fill = "white", color = "black"),
    panel.border = element_rect(colour = "black", fill = NA),
    axis.text.x = element_text(size = 8, hjust = .5,face = "bold"), 
    axis.text.y = element_text(size = 10, hjust = .5), 
    plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
    axis.title.x = element_text(size = 8,face = "bold"),
    strip.text = element_text(size = 10,face = "bold")
  ) +
  geom_text(aes(label = percent(correct_rate, accuracy = 1)), 
            vjust = -0.5, size = 5, color = "black")

ggsave("../paper/p2/p1.pdf", 
       plot = p1,
       device = "tiff",
       width = 17.4,  
       height = 12,    
       units = "cm",
       dpi = 1200,    
       compression = "lzw",  
       bg = "white")


data_raw_long$news_correct_dummy <-  as.factor(data_raw_long$news_correct)



#### Start of Su Code###############################################


data_raw_long$edu3 <- as.numeric(as.factor(data_raw_long$edu_cate))
data_raw_long$exp3 <- ifelse(data_raw_long$exposure %in% 0, 1, 
        ifelse(data_raw_long$exposure %in% c(1,2), 2, 
        ifelse(data_raw_long$exposure %in% c(3:12), 3, NA)))
data_raw_long$over3 <- ifelse(data_raw_long$broad == 0, 1, 
        ifelse(data_raw_long$broad == 1, 2, 
        ifelse(data_raw_long$broad %in% c(2,3), 3, NA)))
data_raw_long$type <- as.numeric(as.factor(data_raw_long$news_category))


library(rstan)
rstan_options(auto_write = TRUE)
rstan_options(threads_per_chain = 1)
options(mc.cores = parallel::detectCores()-1)


# 2PL IRT model

stanCode <- "
data {
  int<lower=1> P;                // Number of respondents
  int<lower=1> J;                // Number of question types (items)
  int<lower=1> N;                // Total number of responses
  int<lower=0,upper=1> y[N];     // Binary response
  int<lower=1,upper=J> type[N];  // Item type for each response
  int<lower=1,upper=P> id[N];    // Respondent ID
  // Covariates (assumed all standardized in R!)
  real treat[N];
  real age[N];
  int<lower=0> ccpmem[N];
  int<lower=0> female[N];
  int<lower=0> college[N];
  real polknow[N];
  real income[N];
  real oversea[N];
  int<lower=0> urban[N];
  real soclevel[N];
  real nationalism[N];
  real ctrust[N];
}
parameters {
  // Latent respondent abilities
  vector[P] theta_raw;
  // Item parameters
  vector[J] alpha_raw;
  vector[J] delta_raw;
  // Covariate effects (betas for each question type or global)
  vector[J] b1;
  vector[J] b2;
  vector[J] b3;
  real b4;
  real b5;
  real b6;
  real b7;
  real b8;
  real b9;
  real b10;
  real b11;
  real b12;
}
transformed parameters {
  // Center and scale for identification
  vector[P] theta = (theta_raw - mean(theta_raw)) / sd(theta_raw);                 // respondent ability
  vector[J] alpha = alpha_raw - mean(alpha_raw);                                   // item difficulty, mean-constrained
  vector<lower=0>[J] delta = exp(delta_raw);                                       // discrimination, positive
  vector[N] eta;
  for (i in 1:N) {
    real linpred = b1[type[i]] * treat[i]
                 + b2[type[i]] * college[i]
                 + b3[type[i]] * oversea[i]
                 + b4 * polknow[i]
                 + b5 * nationalism[i]
                 + b6 * ctrust[i]
                 + b7 * soclevel[i]
                 + b8 * age[i]
                 + b9 * female[i]
                 + b10 * ccpmem[i]
                 + b11 * income[i]
                 + b12 * urban[i];

    eta[i] = delta[type[i]] * (theta[id[i]] - alpha[type[i]]) + linpred;
  }
}
model {
  // Priors: tighter, light-tailed (identification + speed - avoid cauchy)
  theta_raw ~ std_normal();      // global mean=0 sd=1 by construction
  alpha_raw ~ std_normal();
  delta_raw ~ normal(0, 0.5);   // discriminations: lognormal(0,0.5)
  // Regression coefficients: light regularization
  b1 ~ normal(0, 1);
  b2 ~ normal(0, 1);
  b3 ~ normal(0, 1);
  b4 ~ normal(0, 1);
  b5 ~ normal(0, 1);
  b6 ~ normal(0, 1);
  b7 ~ normal(0, 1);
  b8 ~ normal(0, 1);
  b9 ~ normal(0, 1);
  b10 ~ normal(0, 1);
  b11 ~ normal(0, 1);
  b12 ~ normal(0, 1);
  // Likelihood
  y ~ bernoulli_logit(eta);
} 

"

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())


#datalist <- with(data_raw_long, list("y" = as.numeric(news_correct_dummy) - 1,  "type"= type, 
#    "id" = as.numeric(as.factor(id)), "treat" = as.numeric(exp3)-2, "female" = as.numeric(gender),
#    "age" = as.numeric(age), "ccpmem" = as.numeric(partymember),
#    "polknow" = as.numeric(know_pol),
#    "oversea" =  as.numeric(over3)-2,
#    "college" = as.numeric(edu3==4),
#    "income" = as.numeric(incomeFam)-3,
#    "rural" = 2-as.numeric(hukou),
#    "nationalism" = arm:::rescale(nationlism),
#    "soclevel" = arm:::rescale(socialLevel),
#    "ctrust" = arm:::rescale(trust_cen),
#    "N" = length(id), "P" = length(unique(id)), "J" = 4))
    
datalist <- with(data_raw_long, list("y" = as.numeric(news_correct_dummy) - 1,  "type"= type, 
    "id" = as.numeric(as.factor(id)), "treat" = arm:::rescale(as.numeric(exp3)), "female" = as.numeric(gender),
    "age" = arm:::rescale(as.numeric(age)), "ccpmem" = as.numeric(partymember),
    "polknow" = arm:::rescale(as.numeric(know_pol)),
    "oversea" =  arm:::rescale(as.numeric(over3)),
    "college" = as.numeric(edu3==4),
    "income" = arm:::rescale(as.numeric(incomeFam)),
    "urban" = 2-as.numeric(hukou),
    "nationalism" = arm:::rescale(nationlism),
    "soclevel" = arm:::rescale(socialLevel),
    "ctrust" = arm:::rescale(trust_cen),
    "N" = length(id), "P" = length(unique(id)), "J" = 4))

    
inits <- function(){
    list(delta_raw = abs(rnorm(4, 1,.01)))
}

BM01 <- stan(model_code = stanCode, data=datalist, chains=1, iter=10)
BM01 <- stan(fit=BM01, data=datalist, chains=10, iter=2000, cores=10)
BM01 <- stan(fit=BM01, data=datalist, init=inits, chains=8, iter=2500, cores=8, thin=1, warmup=500)

delta <- extract(BM01, pars="delta")$delta
theta <- extract(BM01, pars="theta")$theta
alpha <- extract(BM01, pars="alpha")$alpha
b1 <- extract(BM01, pars="b1")$b1
b2 <- extract(BM01, pars="b2")$b2
b3 <- extract(BM01, pars="b3")$b3
b4 <- extract(BM01, pars="b4")$b4
b5 <- extract(BM01, pars="b5")$b5
b6 <- extract(BM01, pars="b6")$b6
b7 <- extract(BM01, pars="b7")$b7
b8 <- extract(BM01, pars="b8")$b8
b9 <- extract(BM01, pars="b9")$b9
b10 <- extract(BM01, pars="b10")$b10
b11 <- extract(BM01, pars="b11")$b11
b12 <- extract(BM01, pars="b12")$b12



exps <- with(data_raw_long, tapply(as.numeric(exp3), id, mean))
edus <- with(data_raw_long, tapply(as.numeric(edu==4), id, mean))
overs <- with(data_raw_long, tapply(as.numeric(as.numeric(over3)), id, mean))
ages <- with(data_raw_long, tapply(as.numeric(as.numeric(age)), id, mean))

latentb <- apply(theta, 2, mean)
LM01 <- lm(latentb ~ edus + exps + overs + ages)


betas <- as.matrix(coef(LM01))
cefs <- coef(LM01)
pdf("latent.pdf", width=5, height=4)
par(mar=c(3,3,1,1), mgp=c(1.5,0.2,0), tcl=-0.2)
plot(x=ages, y=latentb,type="n", xlab="age", ylab=expression("Individual latent ability ("*theta[p]*")"))
points(x=jitter(ages[edus==1]), y=latentb[edus==1], pch=1, cex=0.3,col="gray")
points(x=jitter(ages[edus==0]), y=latentb[edus==0], pch=19, cex=0.3,col="gray")
curve(cefs[1] + cefs[2]*1 + cefs[3]*3 + cefs[4]*1 + cefs[5]*x, lty=1, add=TRUE, lwd=2)
curve(cefs[1] + cefs[2]*0 + cefs[3]*3 + cefs[4]*1 + cefs[5]*x, lty=2, add=TRUE, lwd=2)
text(x=50, y=cefs[1] + cefs[2]*1 + cefs[3]*3 + cefs[4]*1 + cefs[5]*50, labels="Educ = college degree", adj=1)
text(x=20, y=cefs[1] + cefs[2]*0 + cefs[3]*3 + cefs[4]*1 + cefs[5]*20, labels="Educ < college degree", adj=0)
text(x=55, y=-2, labels=expression(R^2*" = 0.14"), adj=0, font=3)
text(x=55, y=-2.4, labels="N = 2234", adj=0)
dev.off()






quant <- function(x, p=0.025){
    out <- quantile(x, prob=p)
    return(out)
}

namelabs <- c("", expression(alpha["Praising (True)"]), expression(alpha["Praising (False)"]), 
    expression(alpha["Disparaging (True)"]),expression(alpha["Disparaging (False)"]), 
    "", expression(delta["Praising (True)"]), expression(delta["Praising (False)"]), 
    expression(delta["Disparaging (True)"]),expression(delta["Disparaging (False)"]),  "",
    "College degree", "Freq. exposure to political news", "Overseas experience", "Political knowledge", 
    "Nationalism", "Trust in central government", "Perception of social status", "Age", "Female", "CCP member", 
   "Household income", "Urban residency"
   )

    
pts <- c(NA, apply(alpha, 2, mean),  NA, apply(delta, 2, mean), NA, apply(b2, 2, mean), apply(b1, 2, mean), apply(b3, 2, mean), mean(b4), mean(b5),
    mean(b6), mean(b7), mean(b8), mean(b9), mean(b10), mean(b11), mean(b12))
cilo <- c(NA, apply(alpha, 2, quant),  NA, apply(delta, 2, quant), NA, 
    apply(b2, 2, quant), apply(b1, 2, quant), apply(b3, 2, quant), 
    quant(b4), quant(b5), quant(b6), quant(b7), quant(b8), quant(b9), 
    quant(b10), quant(b11), quant(b12))
cihi <- c(NA, apply(alpha, 2, quant, p=0.975),  NA, apply(delta, 2, quant, p=0.975), NA, 
    apply(b2, 2, quant, p=0.975), apply(b1, 2, quant, p=0.975), apply(b3, 2, quant, p=0.975),
    quant(b4, p=0.975),
    quant(b5, p=0.975), quant(b6, p=0.975), quant(b7, p=0.975), quant(b8, p=0.975), 
    quant(b9, p=0.975), quant(b10, p=0.975), quant(b11, p=0.975), quant(b12, p=0.975))


pdf("coefplotaux.pdf", width=6, height=3.5)
par(mar=c(2,5,3,1), mgp=c(1.5,0.2,0), tcl=-0.2, oma=c(0,6,1,0))
plot(0,0,type="n", ylim=c(10,1), xlim=c(-0.7,1.7), xlab="", ylab="", main="", axes=FALSE, frame.plot=TRUE)
points(x=pts[1:10], y=1:10, pch=19, cex=0.6)
segments(x0=cilo[1:10], x1=cihi[1:10], y0=1:10, y1=1:10)
abline(v=0, lty=2)
axis(3)
text(x=-2.1, y=1:10, labels=namelabs[1:10], xpd=NA, adj=0)
text(x=-2.22, y=c(1,6), labels=c("Difficulty parameters","Discrimination parameters"), adj=0, xpd=NA, font=3)
#text(x=-5.55, y=c(1,6,11), labels=c("Difficulty parameters","Discrimination parameters","Covariates"), adj=0, xpd=NA, font=3)
mtext("IRT model parameters", side=3, font=2,cex=1.2, outer=TRUE, line=-1)
dev.off()


pdf("coefmain.pdf", width=6.5, height=4)
par(mfrow=c(2,2), mar=c(1,1,2,0.5), mgp=c(1.5,0.2,0), tcl=-0.2, oma=c(1,12,1,1))
plot(0,0,type="n", ylim=c(3.5,0.5), xlim=c(-1,1), xlab="", ylab="", main="Praising (True)", axes=FALSE, frame.plot=TRUE)
segments(x0=cilo[c(12,16,20)], x1=cihi[c(12,16,20)], y0=(1:3), y1=(1:3))
points(x=pts[c(12,16,20)], y=(1:3), pch=19, cex=0.6)
axis(1, labels=FALSE)
abline(v=0, lty=2)
text(x=-3.25, y=1:3, labels=c("College degree", "Freq. expose to political news", "Overseas experience"), xpd=NA, adj=0)

plot(0,0,type="n", ylim=c(3.5,0.5), xlim=c(-1,1), xlab="", ylab="", main="Praising (False)", axes=FALSE, frame.plot=TRUE)
segments(x0=cilo[c(12,16,20)+1], x1=cihi[c(12,16,20)+1], y0=(1:3), y1=(1:3))
points(x=pts[c(12,16,20)+1], y=(1:3), pch=19, cex=0.6)
axis(1, labels=FALSE)
abline(v=0, lty=2)

plot(0,0,type="n", ylim=c(3.5,0.5), xlim=c(-1,1), xlab="", ylab="", main="Disparaging (True)", axes=FALSE, frame.plot=TRUE)
segments(x0=cilo[c(12,16,20)+2], x1=cihi[c(12,16,20)+2], y0=(1:3), y1=(1:3))
points(x=pts[c(12,16,20)+2], y=(1:3), pch=19, cex=0.6)
axis(1, xpd=NA)
abline(v=0, lty=2)
text(x=-3.25, y=1:3, labels=c("College degree", "Freq. expose to political news", "Overseas experience"), xpd=NA, adj=0)

plot(0,0,type="n", ylim=c(3.5,0.5), xlim=c(-1,1), xlab="", ylab="", main="Disparaging (False)", axes=FALSE, frame.plot=TRUE)
segments(x0=cilo[c(12,16,20)+3], x1=cihi[c(12,16,20)+3], y0=(1:3), y1=(1:3))
points(x=pts[c(12,16,20)+3], y=(1:3), pch=19, cex=0.6)
axis(1, xpd=NA)
abline(v=0, lty=2)
#mtext("Estimated Effects of Major Predictors in the IRT Model, by Information Type", side=3, font=2,cex=1.2, outer=TRUE, line=1)
dev.off()


pdf("coefplot.pdf", width=8, height=5)
par(mar=c(3,3,3,1), mgp=c(1.5,0.2,0), tcl=-0.2, oma=c(0,8.5,1,0))
plot(0,0,type="n", ylim=c(9.5,0.5), xlim=c(-1.1,1.1), xlab="", ylab="", main="", axes=FALSE, frame.plot=TRUE)
points(x=pts[24:32], y=1:9, pch=19, cex=0.6)
segments(x0=cilo[24:32], x1=cihi[24:32], y0=1:9, y1=1:9)
abline(v=0, lty=2)
axis(3)
text(x=-2.1, y=1:9, labels=namelabs[15:23], xpd=NA, adj=0)
#text(x=-1.9, y=c(1,6), labels=c("Difficulty parameters","Discrimination parameters"), adj=0, xpd=NA, font=3)
#text(x=-12.25, y=1, labels=c("Covariates"), adj=0, xpd=NA, font=3)
mtext("Regression estimates of IRT model by socio-political covariates", side=3, font=2,cex=1.3, outer=TRUE, line=-1)
dev.off()






pdf("marginal.pdf", width=8, height=3)
par(mfrow=c(1,3), mar=c(3,1,1,1), mgp=c(2,0.2,0), tcl=-0.2, oma=c(0,3,0,0))


plot(0,0, type="n", frame.plot=TRUE, axes=FALSE, xlim=c(-1,1), ylim=c(0,1), xaxs="i",
    ylab="", xlab="Trust in Central Government")
axis(2, at=c(0,0.5,1))
axis(1, at=c(-0.8,-0.26667,0.266667,0.8), labels=c("Vvery low","Low","High", "Very high"))
#axis(1, at=c(-0.8,-0.26667,0.266667,0.8), labels=FALSE)
#text(x=c(-0.6,0.6), y=c(0.6, 0.4), labels=c("Praising (True)", "Disparaging (False)"), xpd=NA)
colline <- c(rgb(0,0,0), rgb(1,0,0), rgb(0,200/255,150/255), rgb(0,80/255,200/255))
colsim <- c(rgb(0,0,0, alpha=0.02), rgb(1,0,0, alpha=0.02), rgb(0,200/255,150/255, alpha=0.02), rgb(0,80/255,200/255, alpha=0.02))

for(i in c(1,2)){
    for(s in sample(1:16000, size=100, replace=FALSE)){
        curve(arm:::invlogit(alpha[s,i] + delta[s,i] *( 
            b1[s]*1 +  b2[s]*1 +  b3[s]*1 + b4[s]*1 + b5[s]*x + b6[s]*0 +
            b7[s]*1 +  b8[s]*1 +  b9[s]*1 + b10[s]*0 + b11[s]*0 + b12[s]*1) + mean(rp)), from=-1, to=1, add=TRUE, col=colsim[i])
    }
curve(arm:::invlogit(median(alpha[,i]) + median(delta[,i]) *( 
    median(b1)*1 +
    median(b2)*1 + 
    median(b3)*1 +
    median(b4)*1 + 
    median(b5)*x +
    median(b6)*0 +
    median(b7)*1 + 
    median(b8)*1 +
    median(b9)*1 + 
    median(b10)*0 +
    median(b11)*0 + 
    median(b12)*1) + mean(rp)), from=-1, to=1, add=TRUE, col=colline[i])
}




plot(0,0, type="n", frame.plot=TRUE, axes=FALSE, xlim=c(-1,1), ylim=c(0,1), xaxs="i",
    ylab="", xlab="Expose to Online Political News")
axis(2, at=c(0,0.5,1), labels=FALSE)
axis(1, at=c(-0.8,0, 0.8), labels=c("Rare", "Occassional", "Frequent"), adj=0, xpd=NA)
#text(x=c(-0.6,0.6), y=c(0.7, 0.4), labels=c("Praising (True)", "Disparaging (False)"), xpd=NA)
colline <- c(rgb(0,0,0), rgb(1,0,0), rgb(0,200/255,150/255), rgb(0,80/255,200/255))
colsim <- c(rgb(0,0,0, alpha=0.02), rgb(1,0,0, alpha=0.02), rgb(0,200/255,150/255, alpha=0.02), rgb(0,80/255,200/255, alpha=0.02))

for(i in c(1,2)){
    for(s in sample(1:16000, size=100, replace=FALSE)){
        curve(arm:::invlogit(alpha[s,i] + delta[s,i] *( 
            b1[s]*x +  b2[s]*1 +  b3[s]*1 + b4[s]*1 + b5[s]*1 + b6[s]*0 +
            b7[s]*1 +  b8[s]*1 +  b9[s]*1 + b10[s]*0 + b11[s]*0 + b12[s]*1) + mean(rp)), from=-1, to=1, add=TRUE, col=colsim[i])
    }
curve(arm:::invlogit(median(alpha[,i]) + median(delta[,i]) *( 
    median(b1)*x +
    median(b2)*1 + 
    median(b3)*1 +
    median(b4)*1 + 
    median(b5)*1 +
    median(b6)*0 +
    median(b7)*1 + 
    median(b8)*1 +
    median(b9)*1 + 
    median(b10)*0 +
    median(b11)*0 + 
    median(b12)*1) + mean(rp)), from=-1, to=1, add=TRUE, col=colline[i])
}



plot(0,0, type="n", frame.plot=TRUE, axes=FALSE, xlim=c(-1,1), ylim=c(0,1), xaxs="i",
    ylab="", xlab="Overseas Experiences")
axis(2, at=c(0,0.5,1), labels=FALSE)
axis(1, at=c(-0.8,0, 0.8), labels=c("Rare", "Occassional", "Frequent"), adj=0, xpd=NA)
text(x=c(-0.6,0.5), y=c(0.7, 0.4), labels=c("Praising (True)", "Disparaging (False)"), xpd=NA)
colline <- c(rgb(0,0,0), rgb(1,0,0), rgb(0,200/255,150/255), rgb(0,80/255,200/255))
colsim <- c(rgb(0,0,0, alpha=0.02), rgb(1,0,0, alpha=0.02), rgb(0,200/255,150/255, alpha=0.02), rgb(0,80/255,200/255, alpha=0.02))

for(i in c(1,4)){
    for(s in sample(1:16000, size=100, replace=FALSE)){
        curve(arm:::invlogit(alpha[s,i] + delta[s,i] *( 
            b1[s]*1 +  b2[s]*1 +  b3[s]*x + b4[s]*1 + b5[s]*1 + b6[s]*0 +
            b7[s]*1 +  b8[s]*1 +  b9[s]*1 + b10[s]*0 + b11[s]*0 + b12[s]*1) + mean(rp)), from=-1, to=1, add=TRUE, col=colsim[i])
    }
curve(arm:::invlogit(median(alpha[,i]) + median(delta[,i]) *( 
    median(b1)*1 +
    median(b2)*1 + 
    median(b3)*x +
    median(b4)*1 + 
    median(b5)*1 +
    median(b6)*0 +
    median(b7)*1 + 
    median(b8)*1 +
    median(b9)*1 + 
    median(b10)*0 +
    median(b11)*0 + 
    median(b12)*1) + mean(rp)), from=-1, to=1, add=TRUE, col=colline[i])
}


mtext("Pr(Correct=1)", side=2, outer=TRUE, line=0.5)
dev.off()







# comparison between type 1 and type 2 by college and by exposure

plot(0,0, type="n", frame.plot=TRUE, axes=FALSE, xlim=c(-1,1), ylim=c(0,1), xaxs="i",
    ylab="", xlab="Expose to Online Political News")
axis(2, at=c(0,0.5,1), labels=FALSE)
axis(1, at=c(-0.8,0, 0.8), labels=c("Rare", "Occassional", "Frequent"), adj=0, xpd=NA)
#text(x=c(-0.6,0.6), y=c(0.7, 0.4), labels=c("Praising (True)", "Disparaging (False)"), xpd=NA)
colline <- c(rgb(0,0,0), rgb(1,0,0), rgb(0,200/255,150/255), rgb(0,80/255,200/255))
colsim <- c(rgb(0,0,0, alpha=0.02), rgb(1,0,0, alpha=0.02), rgb(0,200/255,150/255, alpha=0.02), rgb(0,80/255,200/255, alpha=0.02))

for(i in c(1,4)){
#    for(s in sample(1:16000, size=100, replace=FALSE)){
#        curve(arm:::invlogit(alpha[s,i] + delta[s,i] *( 
#            b1[s]*x +  b2[s]*1 +  b3[s]*1 + b4[s]*1 + b5[s]*1 + b6[s]*0 +
#            b7[s]*1 +  b8[s]*1 +  b9[s]*1 + b10[s]*0 + b11[s]*0 + b12[s]*1) + mean(rp)), from=-1, to=1, add=TRUE, col=colsim[i])
#    }
for(m in 1:2){
curve(arm:::invlogit(median(alpha[,i]) + median(delta[,i]) *( 
    median(b1)*x +
    median(b2)*1 + 
    median(b3)*1 +
    median(b4)*1 + 
    median(b5)*1 +
    median(b6)*0 +
    median(b7)*1 + 
    median(b8)*1 +
    median(b9)*c(0,1)[m] + 
    median(b10)*0 +
    median(b11)*0 + 
    median(b12)*1) + mean(rp)), from=-1, to=1, add=TRUE, col=colline[i], lty=m)
}}
dev.off()




sfit <- stan_glmer(news_correct_dummy ~ 1+ edu3*factor(type) + ( 1+ edu3*factor(type)|over3) + ( 1+ edu3*factor(type)|exp3) + (1|edu3), data = data_raw_long, family = binomial)


fe <- fixef(sfit)
re <- ranef(sfit)


par(mfrow=c(3,3), mgp=c(1.5,0.2,0), tcl=-0.2, mar=c(1,1,0.5,0.5), oma=c(4,4,4,0))
for(i in 1:3){
    for(j in 1:3){
        plot(0,0,type="n", xlim=c(1,4), ylim=c(0,1), xlab="", ylab="", axes=FALSE, frame.plot=TRUE)
        #abline(h=0.5, col="gray")
        #abline(h=c(0.9,0.4), lty=2)
        if(i==1&j==1){
            text(x=c(2,2),y=c(0.7,0.4), c("Prasing (REAL)", "Praising (FAKE)"), xpd=NA)
            text(x=c(3.5,3.5),y=c(0.6,0.4), c("Disparaing (REAL)", "Disparaing (FAKE)"), col=2, xpd=NA)
        }
        if(i==3& j %in% c(1,2,3)){
            axis(1, at=c(1,2,3,4,5), labels=c("", "","","",""))
            axis(1, at=1, labels="Primary", tck=FALSE)
            axis(1, at=2, labels="Junior", tck=FALSE)
            axis(1, at=3, labels="High", tck=FALSE)
            axis(1, at=3.75, labels="College", tck=FALSE)
            
        }else{
            axis(1, at=1:5, rep("", 5), xpd=NA)
        }
        if(j==1& i %in% c(1,2,3)){
            axis(2, at=c(0,0.5,1), labels=c("0", "0.5", "1"), xpd=NA)
        }else{
            axis(2, at=c(0,0.5,1), labels=FALSE)
        }
        curve(arm:::invlogit(cbind(1, x, 0, 0, 0, 0, 0, 0) %*% (as.matrix(fe) + t(re[[2]][i,]) + t(re[[3]][j,])) + re[[1]][3,]), add=TRUE, lty=1, col=1)
        curve(arm:::invlogit(cbind(1, x, 1, 0, 0, x*1, 0, 0) %*% (as.matrix(fe) + t(re[[2]][i,]) + t(re[[3]][j,])) + re[[1]][3,]), add=TRUE, lty=2, col=1)
        curve(arm:::invlogit(cbind(1, x, 0, 1, 0, 0, x*1, 0) %*% (as.matrix(fe) + t(re[[2]][i,]) + t(re[[3]][j,])) + re[[1]][3,]), add=TRUE, lty=1, col=2)
        curve(arm:::invlogit(cbind(1, x, 0, 0, 1, 0, 0, x*1) %*% (as.matrix(fe) + t(re[[2]][i,]) + t(re[[3]][j,])) + re[[1]][3,]), add=TRUE, lty=2, col=2)
    }
}
mtext("Online Political News Exposure", side=2, outer=TRUE, line=2.5, font=2)
mtext(c("Rare", "Occassional", "Frequent"), at=c(0.85, 0.5, 0.15), side=2, outer=TRUE, line=1)
mtext("Education level", side=1, outer=TRUE, line=1.5, font=2)
mtext("Overseas experience", side=3, outer=TRUE, line=2, font=2)
mtext(c("Rare", "Occassional", "Frequent"), at=c(0.15, 0.5, 0.85), side=3, outer=TRUE, line=0.5)




par(mfrow=c(3,3), mgp=c(1.5,0.2,0), tcl=-0.2, mar=c(1,1,0.5,0.5), oma=c(4,4,4,0))
for(i in 1:3){
    for(j in 1:3){
        plot(0,0,type="n", xlim=c(1,5), ylim=c(0,1), xlab="", ylab="", axes=FALSE, frame.plot=TRUE)
        #abline(h=0.5, col="gray")
        #abline(h=c(0.9,0.4), lty=2)
        if(i==1&j==1){
            text(x=c(3,3),y=c(0.85,0.35), c("Prasing (REAL)", "Praising (FAKE)"))
            text(x=c(3.5,3.5),y=c(0.55,0.7), c("Disparaing (REAL)", "Disparaing (FAKE)"), col=2)
        }
        if(i==3& j %in% c(1,2,3)){
            axis(1, at=c(1,2,3,4,5), labels=c("", "","","",""))
            axis(1, at=1, labels="Primary", tck=FALSE)
            axis(1, at=3, labels="High School", tck=FALSE)
            axis(1, at=4.75, labels="College", , tck=FALSE)
            
        }else{
            axis(1, at=1:5, rep("", 5), xpd=NA)
        }
        if(j==1& i %in% c(1,2,3)){
            axis(2, at=c(0,0.5,1), labels=c("0", "0.5", "1"), xpd=NA)
        }else{
            axis(2, at=c(0,0.5,1), labels=FALSE)
        }
        #for(mm in 1:2) for(s in 1:100)  curve(arm:::invlogit(res[[1]][s,4,] +  (fes[s,] + res[[2]][s,m,]) %*% t(cbind(1, x)) + (fes[s,] + res[[3]][mm,i,])%*% t(cbind(1, x)) + (fes[s,] + res[[4]][s,j,])%*% t(cbind(1, x))), add=TRUE, lty=1, col="gray90")
        for(m in 1:4){
           curve(arm:::invlogit(cbind(1, x) %*% t(aa[[1]][3,]) + cbind(1, x) %*% t(aa[[2]][m,]) + cbind(1, x) %*% t(aa[[3]][i,]) + cbind(1, x) %*% t(aa[[4]][j,])), add=TRUE, lty=c(1,2,1,2)[m], col=c(1,1,2,2)[m])
        }
    }
}
mtext("Online Political News Exposure", side=2, outer=TRUE, line=2.5, font=2)
mtext(c("rare", "occassional", "frequent"), at=c(0.85, 0.5, 0.15), side=2, outer=TRUE, line=1)
mtext("Education level", side=1, outer=TRUE, line=1.5, font=2)
mtext("Overseas experience", side=3, outer=TRUE, line=2, font=2)
mtext(c("rare", "occassional", "frequent"), at=c(0.15, 0.5, 0.85), side=3, outer=TRUE, line=0.5)


pY <- with(data_raw_long, tapply(as.numeric(news_correct_dummy) - 1, type, mean))


pdf("accuracy.pdf", width=6, height=3)
par(mar=c(3,10,3,2), mgp=c(1.5,0.2,0), tcl=-0.2)
plot(0, 0, xlim=c(0,1), ylim=c(4.5,0.5), type="n", axes=FALSE, frame.plot=TRUE,
    xlab="Percentage of Correctness", ylab="", main="Accuracy in News Judgement", xaxs="i", yaxs="i")
segments(x0=rep(0,4), x1= pY, y0=1:4, y1=1:4, lwd=6, col="gray60")
text(x=-0.41, y=0, labels=c("News category"), xpd=NA, adj=0)
text(x=-0.41, y=1:4, labels=c("Praising (True)", "Praising (False)", "Disparaging (True)", "Disparaging (False)"), xpd=NA, adj=0, cex=0.9)
axis(1, at=seq(0,1, length=5), labels=c("0%", "25%", "50%", "75%", "100%"))
abline(v=0.5, lty=2)
dev.off()
#  end of su code #####################




